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ABSTRACT 

The theory and results of modeling for the MSS Control Point Location Error 
Filter are presented. The filter produces the maximum-likelihood estimates for 
average values of the spacecraft position and attitude errors during a single 
scene. The quality of the filter performance is characterized by the maximal 
cross and along- track residual errors for which probability distributions can 
be calculated analytically for a given pattern of control points. The filter 
with an automatic selection of the best set of estimates provides geodetic cor- 
rection at 90% of pixels with residual errors less than 40m for four or more 
control points and the mean-squared measurement errors of the order of 20-25m. 
The same accuracy can be preserved for eight or more control points and measure- 
ment errors of 30-35m. 


INTRODUCTION 

The ground control points (CP) , whose locations are measured on systematically 
corrected imagery and whose true coordinates are known from maps, give highly 
precise information on image displacements at each of the CP’s, The differences 
between true and measured locations provide the input to a filter, which produces 
refined estimates of the spacecraft ephemeris and attitude errors. Then these 
estimates are used for geodetic correction. 

The MSS filter theory, represented i.n Section I, is based on 

1) presentation of image distortions, expressed in Local Space Oblique Mercator 
coordinates, as a linear function of deviations in spacecraft attitudes and 
position (Ref.l), and 

2) recognition of the fact, that MSS processing is limited to a single scene with 

no more than 20 CP’s. It is unlikely that any filter can assess the true time 
dependence in the deviations during a single scene. But we still believe that 
in some cases the MSS filter will be able to produce an reliable estimate of 
average rates. I, 
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An accuracy of these estimates is discussed in Section 2. The covariance 
analysis of the estimate errors shows that image distortions caused by the roll 
and cross-track deviations are so similar that their origins can be determined 
only by near perfect measurements. Thus, the filter is unable to produce an 
reliable estimate of both deviations. At the same time, the filter can provide an 
equivalent estimate for either variable, say, roll, which compensates distortions 
due to both sources. Analogously, for the pitch and along-track deviations. 

The analysis of covariances shows also, that initial uncertainties in rates may 
be reduced only for the equivalent roll + cross- track and pitch + along-track 
rates if there are more than 15 CP’s and the mean-squared measurement errors are 
of the order of 10-15m. So, in cases of few CP’s, that is of interest to us, 
only four estimates should be taken into consideration: for the yaw, radial and 

equivalent pitch and roll deviations. 

Section 3 introduces three global characteristics of filter performance: the 
maximal cross and along-track residual errors, together with combined error in 
distance. These characteristics can be obtained analytically and they establish 
upper levels of errors for any given configuration of CP’s. The final formulae 
for probability distributions are presented^ more details may be found in Ref, 2. 

It is known that pattern, which CP’s form on imagery, have a strong effect on 
filter performance. Examples, given in Section 4, show that one of the most im- 
portant simple characteristics of CP’s distribution is the maximal cross-track 
separation, which has been defined as the maximum of the cross-track distances 
between every pair of CP’s. 

The examples demonstrate also, that for every pattern of CP’s, measurement 
errors, and initial uncertainties in deviations, there is an optimal set of 
estimates, minimizing the residual errors. An approximate algorithm, providing 
the automatic selection of such a set, is described in Section 5. Results of 
modeling indicate that the MSS filter with the automatic selection provides the 
90% average errors less than .5 pixel (40m) for 4 or more CP’s and the mean- 
squared measurement errors of the order of 20-25m, or for 8 or more CP's and 
measurement errors of 30-35m. 


I. THE MSS FILTER EQUATIONS 

Ref. shows that the local SOM coordinates of a frame point, X 

may be . represented as 

X = X + y6, 

o 

where X^ is true coordinates of the point, 5 = (6^, ^2**‘*’ ^ vector of the 

spacecraft position and attitude deviations, and y is a_(2xZ) matrix of the partial 
derivatives (PD), computed at the same point. Now let Z = (Z^,Z 2 ) be the coordinates 
of a CP obtained from a map. We will assume that 

Z = X^ + f , (2) 

where f =(5j^, is a vector of Gaussian measurement errors with zero expec- 

tation values and the covariance matrix R. 


= (X^,X2), 
( 1 ) 
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From (1) and (2), it follows that the measured at a CP displacement, 

AX = X - “Z = - r , (3) 

is also normally distributed with 

E(AX) = E(uT) - E(?) = ]xS 
cov(AX) = E(H^) = R 

Thus, the conditional probability density for AX can be written as 

P = (AX/d) = const'exp j^-J^CAX + y5)^R ^(AX + y6)J 

Let us assume that 6 is constant during a scene. Then a joint density of 
displacements at N control point is 

N 

— -1 — 2 — N n -If 

P(AX ,AX ,...,AX /6) = const •// exp(-% (AX - 

k=l 

r"^(AX^ - y^^)) , 

where upper index k indicates AX and y associated with the k-th CP. 

— A. 

It is known, that the maximum likelihood estimate of 5 (we will denote it as 6) 
is a solution of equation 

VL (T) = 0 , (4) 

where N 

L(d) = In P = 2_ (Ax’^ - r"^(a 5E**^ - y’^i) (5) 

k=l 

and differential operator V is defined in Appendix. 

It is known also that 

E(d) = J (6) 

and in our case (Gaussian conditional density) the covariance matrix of 6 is 

T 

cov(5) = E((6-6)(d- <S) ) = - (vl(6)v^) ^ (7) 

Note from here the summation index k will be omitted. Eq. (4) and (5) yield 

VL = 21 -yd)'^ (M - y?) = 

= *5 r"^(^ - ) = 0 (8) 
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and thus, the solution 6 can be written as 

-1 

6 = M Y , 


where 


Z T' -1 

y R y 


Z - T -1 — 
y R AX 


Noting that 


T V T -1 

~ Ay R 


y R y - - M 


we also have from (7 ) that 


cov(6) = M 


2 ^ 2 


let ^2 be independent with the dispersions ^2 * that case 


and (9) and (10) yield 


^ -1 

6 = M Y 

2—1 

cov(<S) - M 


where M, Y are matrices with elements 

"id ■ "id ■" A- "21 "2d > 


'i " ^ ^"li ^^1 "" 2 "2i ^^2^ 

ft 


Elements of the matrix M ^ will be denoted as m^^ , i.e. M.^ = (m^^) . 


Once 6 is determined, it can be used fo_r geodetic correction. With geodetically 
corrected coordinates of a point being X + p6 , the residual error at the point, 

e = £ 2 ^’ written as 

7 = X + p's - X = X + A- (X + p7) = 
o' 

= y(d-7) (15) 
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Thus, e is normally distributed with zero mean and the covariance matrix 

cov(e) = E jv>(6-6) j" = 

i- _TiT 2 -lT 
= y E {( 6 - 6 )( 6 - 6 r _} V =d^yM y 


( 16 ) 


Eq (16) defines local two-dimensional distribution of the residual ^errors at 
given frame point. It can be used also for detection of outliers , i.e. bad 
measurements at CP's. From (3) it follows, that after geodetic correction, 
measured displacement at the k-th CP's, , can be expressed as 

e. = AX -y6 = y(6-6) +5 
k 

As a sum of two independent Gaussian variables, also is Gaussian with zero 
mean and the covariance matrix. 


cov(e^) =oJyM %’’ + R = Q 


(17) 


The two-dimensional probability density for is represented by the countour 
ellipses 

- —T -1 - A 2 

$(e) = e Q e ” const = A 

It is well-known, that the probability that the 'point' 7^ is Inside the countour 
ellipse is so the k-th CP should be treated as an outlier if 

$(i, ) > , 

K* 2 

where corresponds to a chosen confidence level. For instance, X = 9.21 for 
the 99% confidence level. 

All derived above formulae can be easily generalized to include the case^when 
6 is a slow-changing function of time. Introducing the average deviations and rates 

during a scene, 

a and 3 = (3 » 

we can approximate the deviation at time t as 

+ 


Now the displacement at the k-th CP at time t is 


AX^ = y^( a H- 3t ) + 5 
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and the maximum likelihood function of a, 3 can be written as 

N 

— . — k k: - k- . 

R 

k=l 


L(a, 3 ) = -h ^ (AX^ - y*^(a + gt^)) 


k,— .T 

(AX - y (a - Bt )) . 


From the above, one can obtain that the estimates of a and B (denoted as a and^) 
are given by 


(19) 


A 



a 


Y 

A 

= m7^ 

t 

3 

1 

Y 


where Y is given by (14) and components of Y are 

y'i- p AX )t 

of 

The matrix consists of four submatrices 

* I 


Ml = 


M M 


m' m' ' 


(20) 


where M is defined by (13) and elements of m' and m' ' , m'. . and m', ! , are 

2 IJ 13 

“ij " o2 


f f 

m, . 

IJ 


^2 

= y yo4 i^o-i )t^ 


ll^lj o2 ^21 ^2j 


\ — A — 

We have also that E(a) = a , R(B) = B , and 

A ■ A 2 1 

cov(a,B) = 


( 22 ) 


(23) 


Introducing the estimate of deviation at time t, 6 = a + Bt, we have that 
E(6) = E(a) + t • E(B) = a + Jt = s' 


A 

and the covariance matrix of 6 is 
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cov(^) = E((a-a)(a - a)^) + 2 t E((a-a)(g-i)'^) + 
+ E((6 -e )(t -B = 

« t 2 » » 

= C + 2tC + t C , 

where C, cl and c' * are the submatrices of M. V 

I f i i 



( 24 ) 


(25) 


Further it will be considered that the filter can estimate, at the most, the 
along-track (AT), cross-track (CT) , and radial (RAD) position deviations and 
rates, together with deviations and rates in the pitch (P) , roll (R) and yaw (Y) . 
0^ and will correspond to the cross-track and along-track measurements. 


II. COVARIANCE ANALYSIS 

The covariance matrices, cov(§) and cov (a,^), completely characterize an 
accuracy of estimates, which can be achieved by the filter for a given configuration 
of CP’s and the mean-squared measurement errors and It is well-known, 

that a pattern, which CP’s form on imagery, has a strong effect on filter performance, 
especially for a small number of CP’s. At the same time, our calculations show 
that for N - 10 elements of the covariance matrices insignificantly depend upon a 
distribution of CP’s. For o standard deviation of the estimates are 

approximately proportional to a At the present time, and are not expected 
to be less than 10 and 12m, respectively; the MSS filter will be processing up to 
20 CP’s per scene. 

Table I shows the standard deviations of ^ computed for ct- = 10, = 12m and 

50 (randomly located in a frame) CP’s. Comparison of the standard deviations with 
Initial uncertainties in the spacecraft position and orientation, given in Table 2, 
demonstrates complete inefficiency of the filter in that case. The reason is simple; 
PD with respect to the R and CT deviations, as well as PD with respect to P and AT, 
are almost linearly dependent. As a result, the matrix M is nearly singular, and 
thus Eq. (9) can not give a reliable value of In other words, distortions, caused 
by the R and CT (or P and AT) deviations, are so similar that the difference would be 
revealed only in near perfect imagery by near perfect measurements. 

It prompts not to estimate CT and AP deviations at all, considering the cor- 
responding image dis tori tons as a result of additional fictitious deviations in R 
and P, respectively. Thus, the filter should be treated as a source of appropriate 
geodetic corrections, rather than true estimates. 

The covariance analysis for time-dependent deviations shows that the filter 
is unable to produce reliable estimates of the Y and RAD rates even for N = 50: 

the standard deviations of the estimates are 4-5 times as much as their initial 
uncertainties. 
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At the same time, the filter provides mediocre estimates for combined R + CT 
and P + AT rates when N = 15 - 20 and O- , O 2 of order of 10- 15m. Table 3 
shows such an example for N = 20, = lO, and cr^ “ 12m. Note, that initial 

uncertainties in the R + GT and P 4-'^ AT rates are"^.83 and ,82 yrad/sec (these values 
have been computed by data from Table 2) . 

Despite the fact, that in some favorable conditions the filter can cope with 
these rates, such a case will not be considered below. Being interested chiefly 
in the case of few CP’s we will take into account only estimates of the P,R,Y and 
RAD deviations. 


III. THE MAXIMAL RESIDUAL ERRORS 

The current requirements to geometric correction accuracy are specified in 
terms of .5 pixel 90% of the time. Accordingly, we will evaluate the filter 
performance by the 90% guantile of probability distribution, computed for the 
residual errors which were observed at points of some, say 15x15, grid for randomly 
distributed deviations, measurement errors, and possibly, CP^s locations. Three 
types of the 90% errors may be introduced on two-dimensional grids; for the CT and 
AT components o^f the residuals errors and for the total residual displacements 
e = (e2 + ^ 2 ^^' last characteristic will be referred to as DIST. 

It should be noted, that actually all these characteristics can be obtained 
only by stochastic modeling. At the same time there are two additional global 
characteristics, which can be computed relatively simply: probability distributions 

of the maximal CT and AT errors. These distributions describe errors at the 
worst frame points and thus establish the upper level of possible errors for given 
CP’s. 

A 

Let us introduce the error in the.j-th estimate, A, = 

Eq, (15) may be rewritten as 

"2j (J - 1.2, 3, 4) 


It is known that all PD increase towards the corner points of a frame^ The CT and 
AT errors also reach maximal magnitudes at a corner point j although it is- 

never known beforehand at what specific point. At the corner points only four PD^ 
namely, h21’'^12*^23* *^14’ significant values. Moreover, with an error 

less than .1%, they may be replaced with their maximum values, y (retaining, of 
course, correct sign). Thus, at the corner points 


/V ^ A ^ A 

'l - >'12 *2 + "14 *4 


^ A I ^ A 

^2 ^^21 ^1 ^23 ^3 


Noting that y 2 ^ ^nd PD with respect to Y and RAD, have opposite signs at the 

ends of every scan line, we always can find a corner point, where y^^^ ^£^'^^^14 ^4 
have the same sign (analogously, for W2l^l ^^23 ^3 * indirferent to signs 

of and 02 , we finally have that the maximal (aosolute) CT and AT errors, 

Y2> 
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( 26 ) 


1-1 ' ^ ('‘wM 

'''2 “ I ''21 *1 I I ''23 *3 I 


Here we have omitted the superscript m. 

Fortunately, A., A, and A^, A„ are practically independent and thus these expressions 
may be used separately to'^ derive corresponding distribution and moments (Ref. 2). 


The following are the final formulae for mean value and variance of 
E(Vi> - ^ (Sj + Sp 

Var(Yj^) = (1 - |) (S^ + S^) + 

4- — S. • ( p i»arccos(r) + r - 1) , 

TT i Z 

where 




P = 


2 -1 
^1 ^12 ™22 

"^44 

I 'hi I 





> for i = 1 


J 


(27) 

(28) 


^2 _ 2 2 „-l 

’^21 11 

2 2 ^23 33 


P = 


I 

\/™il ■ "*33 


> for i = 2 


The probability distribution of y. can be written as 


Pr(y^<A) = 


A 2 


1 


J exp (-4^2-) 


2S1 


(^F(B(n + p* t-n * t) + F(B*(n-p • t-n • t)) - ijjdt. 
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where 


n = 1 


B = 


S^*r 


u 


F(u) =-=■ f exp(- t2/2) dt, 
V 2ir V 


and 


^1’ ^2’ given by (29). 


For a given value of error, A, the corresponding probability can be easily 
computed by means of standard subroutines,. Modeling has shown an excellent coincidence 
of theoretical and empirical values of E(y,) and Var(y.). Smirnov's Test also 
demonstrates sufficient coincidence of theoretical and^'empirical distributions. 

Eq. (26) may be used also for evaluation of the maximal residual error in 
distance, which we define as 

■‘n, - <Vi + Tj)'^ 


We could not derive an exact distribution for d . But we have noticed that empirical 
distribution of d^ are similar to Gamma-distribution with the same means and variances 
Because E ( d^) and Var ( can be obtained analytically, we have decided to 


approximate distribution of d^ by Gamma-distribution, which is written here as 


Pr(A) = 


b^r(a) 


rx , 

J "■ 


exp(- ^) 


dll 


where 


a = 


E(d^ )^ 

m 

Var(d^ ) 
m 


(30) 


b = 


Var(d ) 
m 

. E(d^ ) 


(31) 


Because, y ^ and are practically Independent, 


E(d^) = E(yJ) + E(y 2 ) 
Var(d^) = Var(y^) + Var^y^) 


(32) , 

(33) 


In Ref. 2 it is shown that 

2 . _ 2 

r^) 


o 0 0/ 

E(y^) = s, + s,, + ^1^2 ( p afccos (r) + r) 


(34) 
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and 


Var(Y^) = 2(sJ +82)^+ 24 sJs 2 P + 

^ '*’ arccQs(r)) - 

^^tt 2 ( r + p arceos(r))^ (35) 

2 2 

Here and p are given by (29). 2^^* (30)-(35) yield a and b, which are used 

to compute a probability Pr(A) = Pr(A ) for any A by means of a standard subroutine 
for Gamma-distribution . Smirnov’s Test shows sufficient coincidence of the 
approximations and empirical distributions for d ; differences between values of 
errors for corresponding probabilities are less than 5-7%, 


IV. EXAMPLES 

Table 4 presents means, the standard deviations, and the 90% errors for the 
maximal CT and AT errors, together with the 90% errors in distance (BIST). These 
data have been obtained by modeling (M) and analytically (T) for 
and initial uncertainties given in Table 2 (except example 9, where AT = 185m, 

CT = 35m, RAD = 65m, P == R = 120 im'ad, and Y = 35 yrad) . 

From 300 to 600 samples have been used to establish results for each case. The 
examples correspond to five selected configurations of CP’ s , depicted in Fig.l. 
Table 5 describes the examples and shows the mean-squared errors of estimates. 
Examples 1-3 and 7-8 illustrate the fact that for given configuration of CP’s, 
measurement errors, and initial uncertainties, there is an optimal set of estimates 
which provides minimum errors. For distribution A, that set includes P and R for 
distribution C it includes P,R, and Y. Examples 8 and 9 demonstrate also that 
such a set depends upon initial errors in deviations. 

As we already know, the partial derivatives y^^^, y^^, and y^^ represent 

the main effects of the position and attitude deviations on xmage distortions. 

At the same time, there is significant distinction between 1*21 ^14^ ^23* 

when the former are almost constant in a frame, the second increase their 
magnitudes along every scan line. Thus, up to the second order effects, P and R 
estimates do not depend upon position of CP’s in a scene. Roughly speaking, they 
depend upon average CT and AT displacements at all CP’s. On the contrary, to 
detect effects of the Y and RAD deviations, we should observe differences of these 
displacements, so the bigger the CT distances between CP’s, the bigger differences 
in corresponding y^ , and y^Q» the higher an accuracy of the Y and RAD estimates 

Thus, a simple but important characteristic of CP’s distribution is the maximal 
cross- track separation, H, whith we define as the maximum of cross-track distances 
between every pair of CP’s. 

Example 1 shows that for small cross-track separations (H==29.7 km) the Y and 
RAD estimates are absolutely insufficient (457 yrad and 292 m) and, as a result, 
the residual errors are large even for very modest measurement errors. On the 
other hand, even smaller number of CP’s may lead to better results if they are 
’nicely’ separated (example 4 for N = 2 and H = 169 km) . Comparasion of examples 
4 and 5 shows that an along-track shift of CP’s does not affect significantly an 
accuracy of results if cros-track positions are preserved. In addition, example 6 
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suggests that a shift of CP’s as a whole in the cross-track direction towards to 
the frame bounds increases errors. It implies that it is always desirable to 
have CP’s placed symmetrically along the track. 

Analyzing results of modeling, we have noted that the 90% CT and AT errors 
can be approximated as 

Yi = E(y^) + 1.5 yJVar(y^) (36) 

where E(y,) and Var(y^) are given by (27) and (28). We have no explanation of 
that fact^ but it was^verified on a large number of cases which have shown that 
an error of such an approximation usually does not exceed 5%. 


V. AUTOMATIC SELECTION OF ESTIMATES 

As it has been shown, for every pattern of CP’s, measurement errors and 
initial uncertainties, there is an optimal set of estimates which reduces the 
residual errors. Consequently, the filter’s performance can be improved if it 
will automatically select an appropriate set of estimates. Our approximate 
algorithm of selection is based on the fact that y^ and y^ are practically 
Independent, and, bigger maximal errors almost always lead to bigger average errors. 

In our specific case, the a priori known uncertainties in P + AT and R + CT are 
always bigger than errors of corresponding estimates (at least, for mean-squared 
measurement errors less than 40m) . Thus, these estimates always ought to be 
Included in an optimal set. Now, all we need is to compute y^ twice, with and 
without the use of the RAD estimate. In the second case, the standard deviation 
of the RAD estimate must be replaced with the initial mean-squared error. Analogously 
y^ should be computed twice to determine when the Y estimates ought to be employed. 

Table 6 presents the 90% CT, AT, and DIST errors, computed on a 15x15 grid 
as a function of o for various number of CP’s.* Results for each case 

have been established by 300-500 randomly generated sets of CP, measurement errors 
and deviations. CP’s were generated so that the distances between every pair of 
them were not less than 75 km for N ^ 4, 50 km for N = 5,6 and 25 km for N^ 7. 
Additional restriction forbad generation of CP’s on the frame borders. The 90% 

DIST errors also depicted in Fig. 2. 

As one can see, the filter provides geodetic correction with the 90% errors 
less than 40m if cf- — 20m and N — 4. For = 30 m only 8 or more CP’s can 

guarantee that accuracy. Note these results do not include errors due to neglected 
uncertainties in rates. These additional errors, accumulated during 15 sec (i.e. 
with respect to the frame center) can be evaluated as 8.6 m (a) in either direction. 
Being relatively small, they do not affect significantly the total errors. 

It should be pointed out, that the automatic selection only slightly reduces 
the average residual errors. At the same time, it essentially moderates errors 
in relatively rare cases of extremely bad distributions of CP’s. 


* Note, that the CT and AT errors are Gaussian and thus may be described also by 
the corresponding standard deviations. 
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APPENDIX 


a - £ X 1 matrix (vector) 

-T 

a = (a^^ ,a^, . . . ,a^) - transposed vector a 
E(a) - mathematical expectation of a 

cov(a) = E((a - E(a))(a - E(a))^) - covariance matrix of a 

Var(b) - variance of b. 


V = (■ 


95 , 


96, 


• • 96 


-) symbolic differential 


I 


operator defined for U = (u^ , 0 ^ , . . . ,u„ ) as 


T 

7U = 


96 , 


au^ 

96. 


96, 


3u 
A 


If A is a A X it matrix. 



Table 1 


Mean-squared errors in estimates 


(N = 50, = 10, = 12) 



R 

Y 

AT 

CT 

RAD 

yrad 

yrad 

yrad 

m 

m 

m 

2591 

34 

375 

2059 

299 

20 1 
, i 


Table 2 


Initial uncertainties in deviations and rates (la) 



P,R,Y 

AT 

CT 

1 

RAD 1 

deviations 

350 yrad 

550m 

110m 

37m 

rates 

.81 yrad/sec 

. 16m/sec \ 

•065m/sec 

. 65m/sec 


Table 3 


Mean-squared errors in estimates. 



(N = 

20, = 

lOm, 02 = 

12m) 


— 

deviations 


rates 

R + CT 

P + AT 

Y 

RAD 

R + CT 

P + AT 

yrad 

yrad 

yrad 

m 

yrad/sec 

yrad/ sec 

i 3.3 

. , 

4.0 

52.2 

31.1 

.5 

.6 

[.en'.i.' .w.-m* -.a 

. 








Table 4 


The 90% maximal CT, AT, and BIST errors (in meters) 


error 
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Table 5 


Description of Examples in Table 4. 






Table 6 


The 90% CT, AT, and DIST errors (in meters) as a function of N and 
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Figure 2. 

90% errors in distance as a function of measurement errors (a 2 
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